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A finite unbound system which is equilibrium in one reference frame is in general nonequilibrium 
in another frame. This is a consequence of the relative character of the time synchronization in the 
relativistic physics. This puzzle was a prime motivation of the Cooper-Frye approach to the freeze- 
out in relativistic hydrodynamics. Solution of the puzzle reveals that the Cooper-Frye recipe is far 
not a unique phenomenological method that meets requirements of energy-momentum conservation. 
Alternative freeze-out recipes are considered and discussed. 
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I. INTRODUCTION 

Hydrodynamics is now a conventional approach to sim- 
ulations of heavy- ion collisions. Even review papers 
P, 0, 0, 0, S IB] do not comprise a complete list of numer- 
ous applications of this approach. The hydrodynamics is 
applicable to description of hot and dense stage of nuclear 
matter, when the mean free path is well shorter than the 
size of the system. However, as expansion proceeds, the 
system gets dilute, the mean free path becomes compa- 
rable to the system size, and hence the hydrodynamic 
calculation should be stopped at some instant. All hy- 
drodynamic calculations are terminated by a freeze-out 
procedure, while these freeze-out prescriptions are some- 
what different in different models. Moreover, the freeze- 
out prescriptions include recipes to calculate spectra of 
produced particles which are of prime experimental in- 
terest. 

Historically the first method for freeze-out was sug- 
gested by Milekhin Q in the context of the Landau hy- 
drodynamic model of multiple production of particles in 
high-energy hadron collisions [g|. Later, Milekhin's ap- 
proach was criticized by Cooper and Frye 0|. Cooper 
and Frye pointed out that Milekhin's approach does not 
conserve energy and proposed their own recipe of the 
freeze-out. In this paper we would like to discuss a puzzle 
which was in fact a prime motivation of the Cooper-Frye 
approach [§] to the freeze-out in the relativistic hydro- 
dynamics. This puzzle is closely related to the definition 
of the relativistically invariant distribution function as it 
was for the first time advanced by S.T. Belyaev and G.I. 
Budker fijl. 



II. THE PUZZLE 

Let us consider a droplet of matter (for simplicity con- 
sisting of only nucleons) , which is characterized by a total 
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baryon number TV, a total energy E and a total momen- 
tum P, and occupies a volume V. To be precise, we 
assume that this droplet is a closed system. 

Let this droplet be described by an equilibrium distri- 
bution (in configuration and momentum space) 



f(x,p) 



1 



(27r) 3 exp{(p M ^- At ) /T} + 1 



(1) 



in the reference frame charaterized by 4- velocity . Let 
us call this frame as a computation one 1 . This distribu- 
tion is defined in terms of degeneracy of the nucleon g, 
chemical potential /i, temperature T and already men- 
tioned 4- velocity u^. The 4- velocity u M is commonly in- 
terpreted as a velocity with which the droplet moves as 
a whole. We asssume that this distribution is homoge- 
neous in the volume V. The last requirement is an im- 
portant condition of the equilibrium. Therefore, the x 
dependence is in fact absent in Eq. ([I]). In particular, 
distribution function {Ip defines the way how it changes 
under the Lorentz transformation. 

In terms of this distribution function, the conserved 
quantities of the droplet can be expressed as follows. 
First we calculate baryon density (p) and elements of 
the energy- momentum tensor (T^ u ) in the computation 
frame 



P = 

rpfJbV 



d 3 P 



(2) 



PV M = (e + P)^u v - g^P, (3) 



where e and P are the proper energy density and pres- 
sure, respectively. Then we multiply these quantities by 
the volume V and thus obtain 

N = pV, (4) 
E = T 00 V = [(e + P)u°u° - P]V, (5) 

pl = T Oiy = ( £ + p) u O ul y ( 6 ) 

Now we are able to formulate the puzzle. We know 
that (E,P) is a 4-vector, at least this is stated in all 
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textbooks. To be precise, the fact that 

pn = J dV T ^ (7) 

is indeed a 4-vector and that P M is independent of the 
frame 2 (up to a Lorentz transformation), where it is cal- 
culated, is proved, e.g., in Ref. [lH 3 . Then the relation 

P l /E = u*/u° = v l (8) 

should take place, if v l is the velocity of motion of this 
droplet as a whole. As we see from Eqs. ([5]) and ([6]), this 
is not the case. Than the questions arise: what is the 
meaning of the ^velocity and what is the meaning of 
the proper energy density e and the pressure P ? 

Moreover, if we believe that the 4- velocity u M is the ve- 
locity of motion of this droplet as a whole and e is the en- 
ergy density in the droplet-rest frame, we can first calcu- 
late P M in the droplet-rest frame (where w M = (1, 0, 0, 0)) 
and then boost it into the computation frame. Then we 
arrive at another surprising result 

E = eu°V*, (9) 

p i = £u'r, (io) 

where V* is the volume in the droplet-rest frame. Now 
the above puzzle reads as follows. There exists no V* 
which makes Eqs. J9]) and (fTU| compatible with Eqs. ([5} 
and ([6]). This again makes us doubtful about interpreta- 
tion of e and P quantities. 

In fact, precisely the contradiction between Eqs. $9§- 
(fl0|l . on the one hand, and Eqs. {SI-©, on the other 
hand, motivated Cooper and Frye [9J| to suggest their 
recipe for the freeze-out, which just avoids this contra- 
diction rather than resolves it. 



III. RESOLUTION OF THE PUZZLE 

Let us consider what really happens to the equilibrium 
distribution ([I]) under the Lorentz transformation. It is 
convenient to represent this distribution by an ensemble 
of particles as follows 4 

f(x, P ) = ^ 3 (p - Pi(*)) <5 3 ( x - Mm, (ii) 

i 

where x, and Pi are the coordinate and momentum of 
the zth particle, respectively. The Xj coordinates homo- 
geneously populate the volume, V, of the droplet in the 



computation reference frame. By definition of the dis- 
tribution function, all these particles are considered at 
the same time instant t. Integration of this distribution 
function over d 3 p d 3 x with weights 1, p° and p gives 

N = H 1 > E = T,p^ p = Ep- ( 12 ) 

i i i 

respectively, which explicitly demonstrates that (E, P) is 
indeed a 4-vector. 

Let us transform distribution from the computation 
frame JT]), where it is simulated by Eq. ifTTj) . to the rest 
frame of the droplet. To do this, we boost these parti- 
cles with some velocity — v* which certainly differs from 
v l = u l /u° in view of consideration of the previous sect. 
Applying a Lorentz transformation to the ensemble of 
particles (fTTj) . we arrive at 

S 3 (P* - P*i(i*)) <5 3 (x* - x*^*)), (13) 

i 

where quantities marked by * correspond to the rest 
frame of the droplet and are obtaned by the Lorentz 
transformation 5 

t* = xsmht/j + tcoship, 

x* = x coship + tsmhip, y*=y 1 z* = z, (14) 

with tanh i/; = v* . 

We do not call sum lfT3|) a distribution function, since 
all particles are taken at different time instants t*. This 
is a direct consequence of the Lorentz transformation- 
events which are simultaneous in one reference frame are 
not necessarily simultaneous in another one. 

In order to obtain a distribution function from ensem- 
ble i|13p . we should reduce this ensemble to a common 
time, e.g., 

t* = j2n/N, (15) 

i 

where N is the number of particles in this ensemble, cf. 
Eq. (fT2l) . To do this, we should move particles forward 
or backward in time, depending on the sign of t* — t*. 
After this reduction the ensemble (fl3|l already simulates 
a distribution function in the droplet-rest frame: 

f(x*,p*) = J2 S3 (P* -P*i(**)) <5 3 (x*-x*, ( (r)). (16) 

i 

Doing this in general case, we have to take into account 
that particles at time t* in the droplet-rest frame have 
exercised additional (or, vise versa, have not exercised all 
those) interactions as compared to those in the compu- 
tation frame at time t. We will avoid these extra com- 
plications assuming that particle do not interact 6 . This 



2 Experts in the freeze-out prefer to call it as independence of the 
3D hyposurface in the Minkowski space. 

3 See also Ref. [IJ, where this proof is accommodated to the 
problem of freeze-out in nuclear collisions. 

4 Such a representation is extensively used in Ref. |Tll |. 



6 For definiteness, we assume that v* is directed along the x axis. 

6 Moreover, if a system is in a bound state, e.g. a cold nucleus, 
these additional/missed interactions restore equilibrium in any 
reference frame. Therefore, in this paper we consider an inher- 
ently unbound state of the system. 
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case is relevant to the problem of freeze-out. In this case IV. PRACTICAL CONSEQUENCES FOR THE 

FREEZE-OUT 



p%(0 = P*,(*n, 

x%(t*) = x%(fj) H 



(17) 

iCfj/pscoK**-*:). (is) 



i.e. the momentum p*j remains the same, but the coor- 
dinate x*i changes. 

Now we are able to analyze the result of the above 
transformation. Let, for the sake of definiteness, the vol- 
ume V be a Lorentz contracted spherical volume [con- 
tracted with gamma factor 7* = (1— v* 2 ) -1 / 2 ]. The coor- 
dinates Xj(t) homogeneously populate this volume. Since 
the Xi(t) ensemble is taken at the same time instant, 
transformed coordinates x*i(t*) homogeneously populate 
the same but Lorentz "uncontracted" volume, V* = j*V. 
Indeed, the linear transformation {Tj| preserves the spa- 
tial homogeneity of this ensemble. 

When we reduce these coordinates to a common time 
t*, see Eq. 1(18)1 . some high- momentum particles (in the 
droplet-rest frame) leave the V* volume, while the most 
part of low-momentum particles remains in this volume. 
Therefore, the Lorentz transformed distribution becomes 
spatially inhomogeneous and thus even nonequilibrium. 
This is purely relativistic effect, associted with relative 
character of the time synchronization in the relativistic 
physics. This effect is closely related to the fact that, 
if even an unbound system was equilibrium at the initial 
time instant, it becomes nonequilibrium at the next time 
instant because of inhomogeneous expansion of the sys- 
tem. In particular, this is the reason why we failed to find 
a volume V* which makes Eqs. ((SJ) and ((10)1 compatible 
with Eqs. J5]) and ([6]). There exists simply no common 
volume V* for all particles in the droplet-rest frame, if it 
is assumed to be homogeneous in the computation frame. 

Nevertheless, the conventional interpretation of quanti- 
ties entering the equilibrium distribution f2P and the way 
of Lorentz transformation prescribed by it are valid, if a 
considered droplet is an open system surrounded by equi- 
librium medium. Let us transform the distribution (jTTJ) 
in the computation frame by boosting it with the velocity 

v = — u/u°. Now let the volume V be a Lorentz con- 
tracted spherical volume [contracted with gamma factor 
7 = (1 — v 2 ) -1 / 2 ]. Then transformed coordinates Xj(tj) 
homogeneously populate a spherical volume, V = 7V. 
However, in view of discussion in the previous sect., the 
total 3-momentum of the droplet in this "tilded" frame is 
still nonzero, P^O. When we reduce these coordinates 
to a common time t, similarly to Eq. 1(18)1 . some particles 
leave the V volume, but at the same time other particles 
come to this volume from the surrounding medium. Af- 
ter this "particle exchange with the medium" the total 3- 
momentum of the droplet, with already changed particle 
content, becomes really zero, and its momentum distribu- 
tion is really described by Eq. {1} with = (1,0,0,0). 



Let us address the question of observable spectrum of 
particles originating from the frozen out droplet of mat- 
ter. Recollect that this droplet is characterized by the 
total baryon number N, and total energy E, momentum 
P, and volume V in the computation reference frame. All 
these quantities are known from solution the hydrody- 
namic equations. Note that thermodynamic quantities, 
i.e. temperature and baryonic chemical potential are not 
directly known from hydrodynamics. 

From the above discussion we see that first we should 
decide in which reference frame this droplet is equilib- 
rium. There are many possibilities to do this choice. 



A. Freeze-out in Computation Frame 

The first natural choice is that the droplet is equilib- 
rium in the computation reference frame. Then we deter- 
mine the chemical potential fx, temperature T, 4-velocity 
u^, and volume V from Eqs. (J4J) — (J6j) and an equation 
of state (EoS). With all parameters of the distribution 
function (J)) being defined, the invariant spectrum of ob- 
servable particles reads as follows 



E 



d 3 p 



Vp° f( P ,x). 



(19) 



comp. frame 



This spectrum obeys conservations of the baryon number 
N, total energy E and momentum P. Note that this 
recipe of the freeze-out differs both from the Cooper- 
Frye one 0| and from Milekhin's one Q. 

A shortcoming of this recipe is that it is closely re- 
lated to the reference frame of computation. In principle, 
we could do computation in a different reference frame. 
Note that an effective freeze-out in kinetic simulations of 
heavy-ion collisions occurs in the same manner, i.e. the 
history of particle collisions is followed in the reference 
frame of computation. 



B. Freeze-out in Local- Rest Frame 

Another natural construction is as follows. Let us start 
as in the previous sect., i.e. transform distribution from 
the computation frame ((TJ, where it is simulated by Eq. 
(fTT)l . to the droplet-rest frame. To do this, we boost 
the system to the velocity — v* = —P/E which certainly 
differs from v % = u l /u° in view of the previous considera- 
tion. Applying a Lorentz transformation to the ensemble 
of particles described by Eq. (fTT)l . we arrive at ensem- 
ble of particles described by Eq. lfl3]). This ensemble 
still does not simulate a distribution function, since all 
particles are taken at different time instants t*. 

Since we consider freeze-out process, we are not in- 
terested in time instants of these frozen-out particles. 
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Therefore, we artificially attribute the same time instant 
[say, that of Eq. ((15) ] to all particles without changing 
their momenta and coordinates. Then we arrive at an 
equilibrium distribution function ((16) but with 

P*i(0 = P\(ti), (20) 
x*i(f) = x%(i*), (21) 

which differ from (flT) - lfT8) only in definition of x*»(i*). 
This distribution takes place in an "uncontracted" volume 

From the practical point of view, we should solve equa- 
tions 



N = p*V*, (22) 

V = i*V\ (23) 

E = e*u*°V*, (24) 

pi = £ * u «y*, (25) 



supplemented by a EoS, in order to determine /x*, tem- 
perature T*, 4- velocity u*, and volume in terms 
of which the invariant spectrum of observable particles 
reads as follows 

(e%P\ =V* (p^) f*(x,p), (26) 

\ " V J Milekhin 

where f*(x,p) is the equilibrium distribution function 
defined in terms of thermodynamic quantities with su- 
perscript *, cf. Eq. (J). This spectrum obeys conser- 
vations of the baryon number N, total energy E and 
momentum P. This method can be called a modi- 
fied Milekhin's freeze-out, since equations of the origi- 
nal Milekhin's method (H)-© certainly differ from (22)- 
([25) . Precisely this method is used in the model of three- 
fluid dynamics [l2l.[l3|. 

An advantage of this recipe is that the choice of the 
reference frame is unique and independent of the frame 
of computation. However, the entropy is not spectacu- 
larly conserved in this method and thereby requests for 
a special consideration. The entropy conservation can be 
taken into account by replacing Eq. (23) by the equa- 
tion of the entropy conservation, S — <j*V*, where a* 
is the entropy density in droplet-rest frame. This way 
the volume V* becomes an independent variable to be 
determined from this set of equations rather than be- 
ing rigidly defined by the Lorentz contraction factor 7*. 
It was found out that spectra calculated with this addi- 
tional requirement of the entropy conservation coincide 
with those based on Eqs. ((22) -(25) within 1%. It im- 
plies that the entropy is fairly good conserved already 
within the modified Milekhin's method defined by Eqs. 
(22)-(26). 

It is important that two above methods of subsects. 
IIV Al and IIV Bl imply that the global freeze-out hyper- 
surface is in general discontinuous. This hypersurface 
is composed of 3-dimensional pieces Act associated with 
weight (Act n^p^), with which this droplet is represented 
in the total sum over all frozen-out droplets. Here n M 



is the normal 4-vector to the piece Act of the hyper- 
surface. In particular, this weight is Vp° in Eq. (19) 
[n M = (1, 0, 0, 0) in the computation frame] or V*(p^u*^) 
in Eq. ((26) [n^ = u*) . An example of such discontinuous 
hypersurface in (1+1) dimensions is presented in Fig. [1] 
(lower panel). 
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Figure 1: (Color online) Freeze-out hypersurface for hydrody- 
namic evolution the ID step-like slab of nuclear matter. The 
upper panel displays the Cooper-Frye choice for the hypersur- 
face. The lower panel schematically illustrates the modified 
Milekhin's prescription, cf. Eq. ((26), for the hypersurface. 
Arrows indicate local 4-velocities on this hypersurface. This 
figure is borrowed from Ref. [12]. 



C. Cooper-Frye Freeze-out 

The Cooper-Frye hypersurface [§] is constructed on the 
condition that this hypersurface is continuous, see Fig. Q] 
(upper panel) . In the Cooper-Frye approach parameters 
of the distribution function, fx, T and u^, are determined 
from Eqs. (j4)-((6). The invariant spectrum of observable 
particles is expressed as follows 

(e%P\ = Aan^f(p,x), (27) 

V a V / Cooper-Frye 

where is the normal 4-vector to the Act pieces of 
the continuous hypersurface. This formula cannot be 
already associated only with choice of a reference frame. 
It can be done, if n^n^ = 1, i.e. if is time- like. 
However, no frame corresponds to n^n M = — 1. Parts of 
the hypersurface with space-like are unavoidable con- 
sequence of continuity of it. Precisely with these parts 
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connected is a problem of the Cooper-Frye method. If 
n^ 1 < 0, occurrin g at space- like n M , the spectrum of 
Eq. ((271 is negative [ll, Ha]- This is a severe problem of 
the method. Note that above discussed recipes lfT9|) and 
((26)1 do not reveal this problem. 



An important option of the above constructions is 
weather the frozen-out matter is removed from the hy- 
drodynamic evolution or not. This removal is associated 
with certain drain terms, Q and R v , in the r.h.s. of hy- 
drodynamic equations 

= Q, (28) 
d^ v = R v , (29) 

where and T^ v are the baryon current and energy- 
momentum tensor, respectively. An example of such 
drain terms is presented in Ref. fl2T |. 

The Cooper-Frye method unambiguously implies that 
the freeze-out does not affect the hydrodynamic evolution 
of the system, i.e. the frozen-out matter is not removed 
from the hydrodynamic phase: Qcf = and Rq F = 0. 
The Cooper-Frye freeze-out, which is applied in the ma- 
jor part of hydrodynamic calculations now, proceeds in 
the following way. The hydro calculation runs absolutely 
unrestricted. The freeze-out hypersurface is determined 
by analyzing the resulting 4-dimensional field of hydro- 
dynamic quantities on the condition of the freeze-out cri- 
terion being met. 

At the same time, the modified Milekhin's method 
([26|) and the freeze-out in the computation frame lfT9|) 
can be used in both regimes. In both cases the energy 
and momentum are conserved. Examples of the modi- 
fied Milekhin's method with and without removal of the 
frozen-out matter from the hydrodynamic evolution are 
presented in Ref. [12] • The removal of the matter indeed 
affects the system evolution. This influence is illustrated 
in Fig. [2j The freeze-out criterion used in this calculation 
stated that the matter is frozen-out when the local en- 
ergy density e gets lower than 0.4 GeV/fm 3 . The e = 0.4 
GeV/fm 3 characteristic curves calculated with and with- 
out freeze-out turn out to be different. Note that the 
value e = 0.4 GeV/fm 3 is achieved right at the surface of 
the system, if the frozen-out matter is removed. At the 
same time the e = 0.7 GeV/fm 3 characteristic curves, 
which lie quite deep inside the system, remain fairly un- 
affected by the freeze-out. 



V. DISCUSSION 

We considered a puzzle which was in fact a prime moti- 
vation of the Cooper-Frye [i| approach to the freeze-out 
in relativistic hydrodynamics. The puzzle consists in the 
fact that naive calculation of the total energy-momentum 
of unbound equilibrium system does not produce a 4- 
vector and, moreover, depends on the reference frame. 




x[fm] 

Figure 2: (Color online) Characteristic curves, corresponding 
to constant values of the energy density e, for hydrodynamic 
evolution the ID step-like slab of the 4 fm width. Initial con- 
ditions for this slab are constructed on the assumption that 
they are formed by the shock-wave mechanism in head-on 
collisions of two ID slabs at E\ a ^ = 10 A GeV. Thus con- 
structed initial state corresponds to the initial energy density 
£o — 3 GeV/fm 3 . Characteristic curves correspond to e = 0.4 
and 0.7 GeV/fm 3 , calculated with (solid lines) and without 
(dashed lines) removal of the frozen-out matter from the hy- 
drodynamic evolution. This figure is borrowed from Ref. |l2l |. 

We argue that a finite unbound system which is equilib- 
rium in one reference frame is in general nonequilibrium 
in another frame. This is a consequence of the relative 
character of the time synchronization in the relativistic 
physics. Thus, naive assumption that this system is equi- 
librium in any reference frame results in this puzzle. So- 
lution of the puzzle reveals that the Cooper-Frye recipe 
is far not a unique phenomenological method that meets 
requirements of energy-momentum conservation. Alter- 
native freeze-out recipes are considered and discussed. 

The above discussion concerned precisely phenomeno- 
logical methods. Recently microscopic treatments of the 
freeze-out process were advanced based on the Bolzmann 
equation [m [l7| and Kadanoff-Baym equations . It 
was found that these microscopic approaches approxi- 
mately justify Cooper-Frye formula (|27f but only on the 
space-like part of the freeze-out hypersurface (i.e. pos- 
sessing a time-like normal vector). Note that on this part 
of the hypersurface the Cooper-Frye method is very close 
to the modified Milekhin's method ([26) (cf. Fig. [TJ as 
well as to the freeze-out in the computation frame (fl9|l . 
The Cooper-Frye formula on the time-like part of the 
freeze-out hypersurface is not reproduced by these treat- 
ments. Precisely on this part the Cooper-Frye formula 
essentially differs from two above mentioned alternative 
methods and also meets the problem of the negative spec- 
trum. 

Two main conclusions have been drawn from these 
microscopic considerations. First, the frozen-out mat- 
ter should be removed from the hydrodynamic evolu- 



tion. This removal is important for the total energy- 
momentum conservation. This conclusion testifies cer- 
tainly not in favor of the standard Cooper-Frye method. 
Another basic conclusion is that sharp freeze-out at some 
3D hypersurface is a rather rough approximation to the 
spectrum formation, because the freeze-out process is 
fairly extended in space and time. It means that the 
particle emission takes place from an extended 4-volume 
rather than from a 3-dimensional hyposurface as it is 
assumed in all above considered phenomenological meth- 
ods. This conclusion is also supported by kinetic sim- 
ulations, see e.g. [lfl| . Therefore, it makes all above 
phenomenological methods questionable. However, the 
numeric implementation of the microscopic methods de- 
veloped in Refs. [H, US] in 3D hydrodynamic simu- 
lations is highly complicated, because it requires integra- 
tion over future evolution of the system for the calcula- 
tion of the particle emission at fixed time instant. The 
implementation performed in Refs. [l3| is not quite con- 
sistent, since it does not take into account the removal of 



the freeze-out with the hydrodynamic evolution. There- 
fore, we still have to use phenomenological methods of 
freeze-out in actual hydrodynamic simulations of heavy- 
ion collisions. The pending problem is to find out which 
of the phenomenological methods most closely simulates 
results of the microscopic methods. 
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